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Since gene regulatory systems contain sometimes only a small number of molecules, these systems 
are not described well by macroscopic rate equations; a master equation approach is needed for such 
cases. We develop an approximation scheme for dealing with the stochasticity of the gene regulatory 
systems. Using an effective interaction concept, original master equations can be reduced to simpler 
master equations, which can be solved analytically. We apply the approximation scheme to self- 
regulating systems with monomer or dimer interactions, and a two-gene system with an exclusive 
switch. The approximation scheme can recover bistability of the exclusive switch adequately. 



I. INTRODUCTION 

Recently, stochastic nature in small systems has at- 
tracted many attentions [lJ-Q - One of the interesting 
examples of the stochasticity is a gene regulatory sys- 
tem; it has been known experimentally that the gene 
regulatory systems show various phenomena caused by 
intrinsic noise [1, [j| . The gene regulatory systems basi- 
cally consists of genes, RNAs, and proteins. The genes 
could sometimes be activated or repressed by regulatory 
proteins known as transcription factors. The number of 
regulatory proteins is sometimes very small, and there 
are large fluctuations. From a theoretical point of view, 
the gene regulatory systems have been studied a lot us- 
ing Monte Carlo simulations (e.g., [fj, lZ(). In addition, 
in order to gain insights into mechanisms or functions 
of the gene regulatory systems, many analytical studies 
have been done |8Hl8j. For example, if we consider a 
self-regulating gene with monomer binding interactions, 
an exact solution has been already known [12]. When 
one considers more complicated systems, some approx- 
imations are needed. Such approximations have also 
been developed; Fokker-Planck or Langevin equation ap- 
proach dHlO], a variational approach [ll|, [l4. IIUl . and 
self-consistent proteomic field approximation}! 3l] ■ 

A gene regulatory system with only two genes and feed- 
back mechanisms has been studied a lot because it plays 
an important role as a genetic switch; two distinct sta- 
ble states emerges, and they could be switched either 
spontaneously or by external signals. In mathematical 
description for the gene regulatory systems, the RNAs 
are sometimes neglected for simplicity, and only genes 
and regulatory proteins are considered. When we con- 
struct a macroscopic rate equation, in which fluctuations 
in protein copy numbers or gene expression states are 
neglected, the analysis for the rate equation tells us the 
following facts: A system with two mutually repressing 
genes shows a bistability, and cooperative binding of reg- 
ulatory proteins is important for making the bistability 
[T9I |20| . Here, the cooperative binding means that com- 
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binations of two or more proteins need to activate or 
repress genes. The macroscopic rate equation gives mul- 
tiple stable solutions, and each solution corresponds to a 
stable state of the gene regulatory systems, which causes 
the bistability. Hence, for the cooperative binding cases, 
it may be enough to use the macroscopic treatments in 
order to investigate qualitative behavior of the bistabil- 
ity. However, other studies have shown that a so-called 
exclusive switch shows a bistability even when the macro- 
scopic rate equations have only one solution [2l| - |23j . Al- 
though the bistability has been confirmed numerically us- 
ing Monte Carlo simulations, no exact or approximated 
analytical treatment has been proposed yet, to the best 
of our knowledge. 

In the present paper, we develop a new approximation 
scheme for gene regulatory systems. In the approxima- 
tion scheme, there is no need to use continuous descrip- 
tion such as Fokker-Planck or Langevin equations, and 
hence the smallness or discrete properties of the system 
are not neglected. The basic idea of the approximation 
is similar to the "self-consistent proteomic field approx- 
imation" developed by Walczak et al. [l3j]. In the self- 
consistent proteomic field approximation, a joint proba- 
bility for all genes is approximated as a product of prob- 
ability distributions for each gene, and then the interac- 
tions between genes and regulatory proteins can be eval- 
uated 'exactly' in this approximation. In [l3j , only toggle 
switches, which consist of two genes, have been studied; 
as denoted in the discussions in , further approxima- 
tion would be needed for self-regulating systems. We here 
extend the concept of [l3j , and develop a more applicable 
approximation scheme; the interactions between genes 
and regulatory proteins are approximated firstly, and an 
effective interaction is introduced. The new approxima- 
tion scheme would be useful to treat more complicated 
cases, such as the exclusive switch. The new approxi- 
mation scheme enables us to give analytical expressions 
for probability distributions of the numbers of proteins, 
without loss of the discreteness property of the system. 
The effective interactions are estimated self-consistcntly. 
We will demonstrate the usefulness of the approximation 
scheme by using self-regulating systems and the exclusive 
switch without cooperative interactions. 

The present paper is constructed as follows. In Sec. II, 
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we give a brief review of a stochastic model for gene reg- 
ulation. In Sec. Ill, self-regulating systems are studied. 
Section MI CI gives one of the important results in the 
present paper, in which our approximation scheme is pro- 
posed. The proposed approximation scheme is applied to 
the exclusive switch in Sec. VI. Section V is concluding 
remarks. 



II. STOCHASTIC MODEL FOR GENE 
REGULATION 

We here briefly review the basic biology of genetic regu- 
latory system and a simplified stochastic model, for read- 
ers' convenience. 

A gene regulatory system consists of many compo- 
nents, such as genes, RNAs, and proteins. The tran- 
scription of a gene is initiated by a binding of RNA poly- 
merase to a promoter site of the gene in the DNA. The 
binding of regulatory proteins (or molecules), so-called 
transcription factors, can sometimes regulate the tran- 
scription initiation. These regulatory proteins bind to 
own target operator sites, and they sometimes act as re- 
pressors (which repress the transcription) or activators 
(which enhance the transcription) of the transcription. 
When the RNA polymerase binds to a gene, the gene 
sequence is copied into a messenger RNA (mRNA), and 
the mRNA is translated into a protein molecules by a 
ribosome enzyme complex. The produced proteins are 
important to determine the phenotypic behavior of the 
cell. In addition, regulation of transcription is one im- 
portant way of controlling the phenotypic behavior, and 
sometimes the produced proteins can become regulatory 
signals for genes. 

Although all of the above reactions would be important 
for the gene regulatory systems, the mRNA is sometimes 
neglected in stochastic modelling for simplicity. That is, 
the translation from mRNAs to proteins are straightfor- 
ward, and then we assume that an activated gene directly 
increases the number of proteins. In addition, we con- 
sider that a repressed gene cannot produce any proteins, 
which makes analytical treatments much simpler [12| . 

In the present paper, all regulatory proteins act as re- 
pressors. If regulatory proteins are not binding to a gene, 
then we call a state of the gene as 'ON' state; if not, the 
gene is in 'OFF' state. A gene in OFF state cannot pro- 
duce any proteins, as we assumed above. 



III. SELF-REGULATING SYSTEM 

In this section, we will explain a new approximation 
scheme using a simplest model, i.e., a self- regulating sys- 
tem. Exact solutions for the self-regulating system with 
monomer interactions have already been known. After 
reviewing the exact solutions, we will propose a new ap- 
proximation scheme. The new approximation scheme will 




FIG. 1: A schematic illustration of the self-regulating gene. 



be applied to the self-regulating systems with monomer 
and dimer interactions, respectively. 



A. Model 

At first, we give a brief explanation for a self- regulating 
system. In the self-regulating system, there is only one 
gene, and it produces proteins. The produced proteins 
are considered as regulatory proteins for the gene, and 
the regulation is a repressed one. In this sense, there 
is a self-regulation mechanism. Figure [Tj shows the self- 
regulating system. When the gene is in ON state, it 
produces proteins with rate g. The degradation rate of 
the regulatory proteins is k. The regulatory proteins can 
bind the gene with rate function H(n), where n is the 
number of 'free' regulatory proteins. The function H(n) 
can be a complicated function of the regulatory proteins; 
e.g., H(n) — hn for monomer interactions, and H(n) — 
hn(n — l)/2 for dimer interactions, where h is a rate 
for the binding. / is the rate with which the regulatory 
protein is released from the repressor site of the gene. 



B. Exact solution for monomer interactions 

We here consider a simplest interaction case, i.e., a 
monomer interaction case. Hence, H(n) in Fig.[Tjis writ- 
ten as hn, as discussed in Sec. MI A1 For the monomer in- 
teraction cases, exact solutions have already been known 
[f| • In order to compare our approximation scheme, 
which will be proposed in Sec. MI CI we here briefly re- 
view the exact solutions. 

In order to make analytical treatments simpler, one 
assumption should be included @; i.e., one of the pro- 
teins in ON state is inert, and then the protein cannot 
be degraded or repress to the gene. Hence, there are 
a little difference between usual stochastic simulations 
and this analytical treatment. However, it has already 
been discussed that this assumption alter only for lower 
numbers of proteins, and actually it gives quantitatively 
good results @. Hence, we here employ this assump- 
tion. Figure [5] shows the transition scheme for the usual 
stochastic simulations and the analytical treatment. a n 
and j3 n correspond to probabilities with which there are n 
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FIG. 2: Transition scheme for simulations and analytical 
calculations for monomer interaction cases. ot n and /3 n cor- 
respond to probabilities with which there are n regulatory 
proteins for ON and OFF states, respectively. In the analyt- 
ical calculations, one of the proteins is considered as an inert 
one when the gene is in ON state, and the inert protein is 
also included in a„; the number of 'free' regulatory proteins 
in ON state is n — 1. 



regulatory proteins for ON and OFF states, respectively. 
In the usual stochastic simulations, the degradation rate 
of the proteins, i.e., the change from a n to a„_i, is pro- 
portional to the number of proteins, n. In contrast, the 
above assumption means that the rate from a n to a n _i 
in the analytical treatment is proportional to n — 1, not 
to n. 

In this assumption, the master equations are given as 



da' 

-jT =?K-l - "'J + k [( n + l)«n+l _ na 'nl 
- hna' n + fp n , 



+ hna' n - //?„ 



(1) 
(2) 



where oi n is a probability with which there is n 'free' 
regulatory proteins for ON state; a n+ \ = a' n . Note that 
the inert protein is not a 'free' regulatory protein, and 
a' n does not include the inert protein. 

In order to solve eqs. ([TJ and ([2]), it is useful to define 
the following generating functions; 



(3) 
(4) 



n=0 



Using the generating functions, various information 
about the self-regulating system can be obtained. For 



example, the number of 'free' regulatory proteins is 
d n 



dz' 



: a'(z) 



(5) 



z=0 



The probability with which the gene is in ON state is 
given as a'(l); the number of total proteins in the system 
is given as 



da'(z) 



dz 



+ 1 x q'(1) 



dj3(z) 



z=l 



dz 



(6) 



z=l 



where the second term in r.h.s means a contribution from 
the inert protein in ON state. 

Eqs. ([1]) and ^ can be rewritten as two differential 
equations in terms of the generating functions; 



da(z) 
dt 

dp(z) 
dt 



(2-1) 



ga (z) - k— 



da' 

hz— + fP(z) 



dz dz 



(7) 
(8) 



After some calculations, 'stationary' solutions for the 
generating functions are obtained as follows [1, [l2[ : 

a' cx (z) =A CX F(l + a cx , 1 + & ox , N cx {z - ^ x )), (9) 



f 



(z-1) 



ga' ex (z) 



da' ex (z) 



dz 



h da' cx (z) 



(10) 



where the super-script 'ex' means 'exact solutions', and 
k 



k + h' 



N c 



h 



a cx = {, &«=_Z_ + (l-*g*)JV~ (11) 
k k + h 

A cx is the normalization constant, which is determined as 
a' cx (l) + /3 GX (1) = 1. F(p, q, r) is the Kummer confluent 
hypergeometric function, 



(P)n 



(q)n n\ ' 



(12) 



where (p) n = p(p + l)(p + 2) • • • (p + n — 1). 

Details of the characteristics of the exact solution are 
written in Q . For example, the probability distributions 
for the numbers of 'free' regulatory proteins are 

/lex / 1 I cx\ 

" ~ n\ ( ' (1 + 

x F(l + a cx + n, 1 + 6 0X + n, -N cx z° x ), (13) 
A cx 
1 



B ex = — 



((1 - z° x )N ex 



(N ex y 

71 



x V ■ ' n F{l + a cx + n, 1 + b cx + n, -N cx z% x ) 



+V 



(1 + b cx Jn 

(N ex ) n (a ex ) n 
n\ (6 ox )„ 



F(a cx + n, b cx + n, -N ex z% x ) 



(14) 
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where f — f /(k + h). 



C. Approximation scheme 



In Sec. IIIIB1 we analyzed the self-regulating system 
with monomer interactions. In the system, the interac- 
tion factor, H(n), is given as hn, and actually this sim- 
ple form of the interaction enables us to obtain the exact 
solutions. If we consider different types of interactions, 
such as dimer ones, exact solutions have not been known 
yet. 

Here, we propose a new approximation scheme. The 
key of the approximation scheme is to use "an effective 
interaction" . Although the new approximation scheme is 
similar to the self-consistent proteomic field approxima- 
tion in 13], the new one is more applicable; it is applica- 
ble even to the self-regulating systems or exclusive switch 
cases, as shown later. The effective interaction means 
that the interaction factor in Fig. [TJ H(n), is replaced 
as a scalar value; i.e., H(n) — h. Here, the effective in- 
teraction h is not a function of the regulatory proteins. 
Hence, the master equation for this approximated system 
is written as follows: 



da' 



k[(n+l)a' n+1 -na' n ] 
= + k[(n+l)P n+1 -n/3 n ] 



- ha' n + f0 n , (15) 

dt 

+ ha' n - f/3 n . (16) 



Because the interaction factor H(n) = h has a simple 
form, the analytic solution can be easily calculated by 
using the generating function approach. Putting the left- 
hand sides of Eqs. (fT5|) and (fTo| as zero and rewriting 
Eqs. (fT5|) and (|16p in terms of the generating functions 
a'(z) and (3(z), stationary solutions for the generating 
functions are obtained as follows: 



a'(z) =AF{a,b,N(z - 1)), 



(17) 



P(z) = ( 1 + j ) AF(a - 1, 6 - 1, N(z - 1)) - a(z), 



where A = f / (/ + h) and 



N = r- 



f_ 

fc' 



1 + 4- , b 



1 



f + h 



(18) 



(19) 



A remaining task is to determine the effective inter- 
action h. For the self-regulating system with monomer 
interactions, the binding of the regulatory proteins oc- 
curs only when the system is in ON state. Hence, the 
number of proteins, which can be attached to the bind- 
ing site, should be equal to the number of free proteins 
for ON state. 



According to the following discussions, we here set the 
effective interaction h as 



h = h{n) a i 



(20) 



where (n) a i is the expectation of the number of free reg- 
ulatory proteins under a condition that the gene is in ON 
state (conditional expectation). Because it is possible to 
evaluate the conditional expectation using the generating 
function (Eq. (TTT])) as 



(n)a' = -^-a'(z) 
a'(l) az 



z=l 



g(k + f) 

Hf+f+hy 



(21) 



we obtain the following self-consistent equation by insert- 
ing Eq. HUD; 



h(n) a ' = h- 



g(k + f) 



Kf + / + %>«') 

Solving Eq. ([22]) . we finally obtain 



(22) 



h(n) a > 



-(k 2 + kf) + ^(k 2 + kf) 2 + 3khg(k + f) 



2k 



(23) 



Once the effective interaction h is determined, all sta- 
tistical properties related to the number of regulatory 
proteins are evaluated from the generating functions 
(Eqs. (|17p and ((HI)). For example, the probability dis- 
tributions for the numbers of free proteins are 



a' =AN n ^F(a + n,b + n, -N), 
x F(a - l+Ti, 6-1 + n, -N) 



(24) 



(25) 



D. Results for monomer interactions 

For monomer interaction cases, the exact solutions are 
obtained. Hence, we here compare the exact results 
and approximate results obtained by the approximation 
scheme. 

For the comparison, we here introduce rescaled param- 
eters as follows 1 61: 



f 
k' 



f 

h' 



ad 



g_ 

2fc' 



(26) 



and for simplicity, we set k = 1 in all numerical evalu- 
ations. These rescaled parameters are helpful to under- 
stand properties of the genetic switch. The parameter 
A a( j characterizes the synthesis/degradation processes, 
and large A a( j would give a large average number of pro- 
teins. X eq is related to the equilibrium constant of the 
binding/unbinding process. Finally, w is a parameter 
called "adiabaticity parameter" . u measures how rapidly 
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FIG. 3: Comparison between results from the exact solu- 
tions and the approximation scheme, (a) The probability dis- 
tributions of the number of proteins. Case A: X eq = 10.0, 
X ad = 10.0, w = 0.01. Case B: X cq = 50.0, X ad = 50.0, 
oj = 10.0. For Case A, it is difficult to see the difference be- 
tween the exact solution and approximate solution, (b) The 
average number of proteins. X eq = 10.0 and X a d = 10.0, and 
only the value of uj was changed. 



the gene can equilibrate in a gene state. If ui is small, the 
synthesis/degradation behaves almost like an indepen- 
dent birth and death process, and there would be two 
peaks corresponding to the binding/unbinding states^re- 
spectively. For details of these parameters, e.g., see 

Firstly, the probability distributions of the number 
of protein were compared. Figure [3][a) shows the re- 
sults. Here, we performed two cases: In Case A, we set 
X cq = 10.0, X ad = 10.0, uj = 0.01; in Case B, X cq = 50.0, 
X ad = 50.0, u> — 10.0. For Case A, the exact solution 
and the approximate one give a good agreement, and 
it is difficult to see the difference. Although there are 
quantitative differences between the exact and approx- 
imate solutions for Case B, the approximation scheme 
gives qualitatively good result despite of the rough ap- 
proximation. Figure [3jb) shows the average number of 
proteins (Eq. (JB])) for various values of u> when X eq = 10.0 
and X ad = 10.0. If uj is small, the approximation scheme 
gives quantitatively good results. Even in the large uj 



Stochastic simulation 




o- 



3fc 



4k 



A) 



Analytic solution 



3k 



4k 



FIG. 4: Transition scheme for simulations and analytical 
calculations for dimer interaction cases. 



case, the difference between the exact and approximate 
results are less than 1. 



E. Results for dimer interactions 

As a second example, a self-regulating system with 
dimer interactions are considered. In this case, the tran- 
sition scheme for analytical calculations are different from 
the monomer interaction cases; see Fig. |4l In this case, 
the master equations are the same as Eqs. ([T5|) and (|16|) . 
but the effective interaction should be set as 



(n(n - l)) a > 



(27) 



and we should interpret a' n as a n +2 = cJ n . 

Using the similar procedure written in Sec. 1111 CI the 
effective interaction h is obtained by solving the following 
self-consistent equation: 



(n(n - l)) a ' 



1 d 2 



2a'(l) dz 



(28) 



z=0 



Since it is a little complicated task to obtain the analyt- 
ical expression for the effective interaction h, we numer- 
ically solved the self-consistent equation (Eq. (|2"g|) b 

For the dimer interaction cases, we have not obtained 
any exact solution. Hence, results of the approximation 
scheme were compared with those of the Monte Carlo 
simulations. 

Figure [S] shows the probability distributions. In Case 
A, we used the following rescaled parameters: X cq 
1 0.O. X ad = 10.0, ui = 0.01; in Case B, X cq = 50.0, 
X ad = 50.0, uj = 10.0. The numbers of the Monte Carlo 
steps are over 10 7 for Case A, and 10 s for Case B. Obvi- 
ously, the approximation scheme gives qualitatively good 
results; although the shapes of the distributions and the 
positions of peaks are slightly different, the number of 
peaks are the same as the Monte Carlo results. In addi- 
tion, the average number of proteins are almost the same 
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FIG. 6: A schematic illustration of the exclusive switch. 



as the Monte Carlo results: for Case A, (n) = 1.2 in the 
Monte Carlo simulation and (n) = 1.3 in the approxima- 
tion scheme; for Case B, (n) = 19.0 in the Monte Carlo 
simulation and (n) = 19.0 in the approximation scheme. 



IV. EXCLUSIVE SWITCH 

Next, we consider a more complicated case, i.e., an 
exclusive switch [21-23]. The exclusive switch consists 
of two genes, gene 1 and gene 2. The two genes have 
overlapping promoter sites, and the binding of one of 
the regulatory proteins prevents the binding of the other 
regulatory proteins. 

Here, the interaction between the binding sites and 
proteins is assumed to be monomer interactions. Because 
the interaction is not cooperative bindings, the macro- 
scopic rate equation gives only one solution (2lT[23l |. Al- 
though there is only one solution, it has been shown that 
the exclusive switch can play as a switch. In the exclu- 
sive switch, stochastic effects make the bistability even 
without cooperativity between the regulatory proteins. 

In order to study the exclusive switch analytically, 
master equations for a joint probability P(n\, U2, si, S2) 
should be constructed; rij G N is the number of proteins 



for gene i, and s, £ {ON, OFF} indicates the gene state. 
In general, the master equations for multiple gene cases 
are very complicated, and it could be difficult to obtain 
numerical solutions if the number of genes is large. 

However, in our approximation scheme, the effective 
interaction is used, and it enables us to reduce the equa- 
tions to be solved. Because we consider only the effec- 
tive interaction, the probability P(ni, n.2, Si, S2) can be 
expressed as P (n\, s\) x P(n%, 82)- For example, the mas- 
ter equation for gene 1 is written as 



da 



'(i) 



dt 



'(i) 



(i) 



dt 



+ fcW[(n+l)/?i 1 ] 1 -n/3W] 



(29) 



(30) 



where the super-script '(1)' indicates gene 1. Master- 
equations for gene 2 can be obtained in the similar way. 
Using the same discussion in the self-regulating systems 
in Sec. Hill the generating functions for genes i 6 {1,2}, 
Qi'W(z) and /?M(z), are derived. 

The effective interaction should be chosen as fol- 
lows. The transition of gene 1 from ON state to OFF 
state can occur only when the gene 2 is in ON state, and 
the effective interaction includes only a contribution 
from the free proteins 2 in ON state. Note that gene 1 
does not know whether gene 2 is in ON state or OFF 
state, different from the self-regulating system discussed 
in Sec. Mil in the self-regulating system, the gene knows 
the own state. Hence, the evaluation of the effective in- 
teraction is slightly different from the self-regulating sys- 
tems. The conditional average of the number of free pro- 
teins 2 is given by ( n^) a , = (da' (2 \z)/dz\ z=1 )/a' (2 \l), 
as discussed in Sec. MI Cl The probability P(2on), with 
which gene 2 is in ON state, is calculated as P(2on) = 
a'^\l). Defining {n^)' a , = {n^) a , P(2 ON ), the effec- 
tive interaction should be written as 



W> =hW(nW)' a ,. 



(31) 



According to the above discussions, we finally obtain the 
following self-consistent equations; 





/l (2) 5 (l) / (l)( fc (l) 




fc(l)(/(l) 


+ hW(n( 2 ))' a ,)(kW 4 








(32) 










h {l) g {2)j{2)( k (2) 


f/ (2)) 



/fc(2)(/(2) +M 2 )(n( 1 ))' Q ,)(fc< 2 ) +/( 2 ) +h,w(nwy a ,y 

(33) 

By solving the above self-consistent equations, we obtain 
M 1 ) and h^ 2 \ We here solve them numerically. 
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FIG. 7: Joint probability distributions for the exclusive switch, (a) Monte Carlo results for X^q = = 25.0, X^ — 

X^ 1 = 25.0, and = u/ 2 ^ = 0.1. (b) Approximate results. All parameters are the same as (a), (c) Monte Carlo results for 
X^ = 40.0, X ( cq = 30.0, X$ = 30.0, X ( *> = 20.0, u; (1) = 10.0 and oj (2) = 20.0. (d) Approximate results. All parameters are 
the same as (c). 



We can immediately obtain the probability distribu- 
tions for each gene using the approximation scheme. In 
order to reconstruct the joint probability distribution 
for genes 1 and 2, we need more calculations as fol- 
lows. Firstly we calculate conditional probabilities for 
the number of free proteins for gene i (i G {1,2}) as 



a'W(l)' 



(34) 
(35) 



Secondly, because of the approximation scheme, the joint 
probability distribution should be evaluated as 

P(n 1 ,n 2 , Ion, 2 n) = P(1on, Zon)®^-!®^-!, (36) 
P(ni,n 2 ,l N,2oFF) = P(1on, 2 OF )a^ ) _J^, (37) 

P(ni,«2,l0FF,20N) = P(l0FF,2 N)/3i 1 1 ) "n 2 2 -l' ( 38 ) 

where P(1on,2on) is the probability with which gene 
1 is in ON state and gene 2 is in ON state, and so on. 
Note that a' means only the number of 'free' proteins; for 
monomer interaction cases, the difference between a and 
a' is only one inert protein. In addition, the probability 



with which both genes 1 and 2 are in OFF state is zero; 
-P(IofF) 2off) = 0, because of the exclusive settings. 

Taking the exclusive settings into account, the 
marginal probabilities are calculated as follows: 



-P(1on, 2on) + P(1ofF) 2on) = P(2on)> 
P(1on, 2off) = -P(2off), 
P (Ioffj 2on) = P(1off), 
k P(loN, 2on) + P(1on, 2off) = P(1on), 



(39) 



and then 



' P(l ON , 2 ON ) = P(1on) - P(2off), 
P(1on,2off) = P(2off), 
P(1off,2on) = -P(Ioff), 
P(1qff, 2qff) = 0. 



(40) 



The marginal probabilities, such as P(1on), can be evalu- 
ated by using the generating function a'W (z) and ( z ). 
Finally, we can construct the joint probability distribu- 
tion as 



P(rai, ra 2 ) =P(ni, m, Ion, 2on) 

+ P(m,n2, 1off,2on) 
+ P(m, n,2, Ion, 2qff)- 



(41) 
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We here note that the probabilities P(1on,2on), calcu- 
lated using the above procedures, may become negative 
for some cases; i.e., -P(Ion) > -P(2off) for some choices 
of parameters g( % ' , few, and /W. In these cases, other 
procedures to estimate the joint probabilities are needed. 
In the following numerical experiments, only the former 
cases (P(1on) < -P(2off)) are treated. 

Figure [7] shows the joint probability distributions. Fig- 
ures[7fa) and (c) are Monte Carlo results, and Figs.[7]Jb) 
and (d) are results of the approximation scheme. As 
in Sec. IIIID1 we used the rescaled parameters, and set 

k (l) = k (2) = L j n Fig ^ [Jj^ and ( b ) ; wc uged the 

parameters X$ = X$ = 25.0, = X { ^ = 25.0, 

and = u {2) = 0.1; for (c) and (d), X$ = 40.0, 
X$ = 30.0, = 30.0, = 20.0, = 10.0 

and u/ 2 ) = 20.0. The numbers of the Monte Carlo steps 
are over 10 s for Fig. [7]Ja), and over 10 9 for Fig. [TJc). 
Although Fig. EJd) does not show the correlated behav- 
ior seen in Fig. [7jc) because correlations between gene 
1 and gene 2 are largely neglected in the approximation 
scheme, one could say that the approximation scheme 
gives qualitatively good results; the characteristics of the 
peak structure can be recovered adequately despite the 
rough approximation. In Fig.[7£b), the bistability due to 
the exclusive settings is recovered well. 

V. CONCLUDING REMARKS 

In the present paper, we developed the approxima- 
tion scheme for gene regulatory systems. We firstly ap- 
plied it to self-regulating systems. The approximation 
scheme gives qualitatively good results; the characteris- 
tics of peak structures can be recovered well. In addition, 
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due to the extension of the basic idea of the effective 
interactions, we can naturally apply the approximation 
scheme even to the exclusive switch, and the bistability 
of the exclusive switch without cooperative interactions 
is successfully recovered. 

In contrast to the Fokker-Planck or Langevin ap- 
proach, the approximation scheme proposed in the 
present paper does not neglect discrete properties of sys- 
tems. In addition, because we can rewritten the joint 
probability for all genes as a product of probability dis- 
tribution for each gene, the dimensions of the problems 
are reduced largely. 

The approximation scheme developed in the present 
paper would be a very crude one; it cannot treat corre- 
lated characteristics between genes. However, approx- 
imate analytic expressions are immediately obtained, 
and qualitatively good results are given despite the 
crude approximation. Since Monte Carlo simulations are 
sometimes time-consuming and need high computational 
costs, it would be beneficial to study such approxima- 
tion scheme in order to obtain qualitative pictures for 
the probability distributions. In addition, developments 
of analytical treatments would be helpful to gain insights 
and understandings for the regulatory systems 
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